Multidimensional Alternating Kernel Method for cortical layer segmentation in 3D reconstructed histology

The neocortex of the brain can be divided into six layers each with a distinct cell composition and connectivity pattern. Recently, sensory deprivation, including congenital deafness, has been shown to alter cortical structure (e.g. the cortical thickness) of the feline auditory cortex with variable and inconsistent results. Thus, understanding these complex changes will require further study of the constituent cortical layers in three-dimensional space. Further progress crucially depends on the use of objective computational techniques that can reliably characterize spatial properties of the complex cortical structure. Here a method for cortical laminar segmentation is derived and applied to the three-dimensional cortical areas reconstructed from a series of histological sections from four feline brains. In this approach, the Alternating Kernel Method was extended to fit a multi-variate Gaussian mixture model to a feature space consisting of both staining intensity and a biologically plausible equivolumetric depth map. This research method• Extends the Alternating Kernel Method to multi-dimensional feature spaces.• Uses it to segment the cortical layers in reconstructed histology volume. Segmentation features include staining intensity and a biologically plausible equivolumetric depth map.• Validates results in auditory cortical areas of feline brains, two with normal hearing and two with congenital deafness.


Background
Much of the cerebral cortex, called neocortex, can be divided into six layers traditionally identified by Roman numerals with increasing depth from I to VI ( Fig. 1 ).Laminar structure and thickness are closely related to information processing in the cerebral cortex, which underlies many cognitive processes including sensory perception and behavior [ 1 , 2 ].However, an exact link between the cortical layer structure and cortical function is not settled yet.For example, differences in the hearing experience are known to alter activation patterns in the auditory cortex [ 3 , 4 ].Prior histological studies have shown that congenital deafness alters auditory cortex laminar thickness with inconsistent results [5][6][7][8][9] .Thus, finding a reliable method for cortical layer segmentation could potentially resolve this problem.
It was recently shown that cortical layers III, V and VI mostly drives cortical thicknesses diversity in sensory, motor and association areas [ 10 , 11 ].This was done through whole-brain laminar segmentation of a silver-stained human brain (BigBrain) histological atlas via a convolutional neural network trained on 51 manually segmented regions.A different approach took advantage of the rich amount of information available in diffusion MRI [12] .Layers were separated through k-means clustering based on 31 diffusivity and orientation parameters.Classification did not use spatial or depth information but instead relied on the high-dimensional parameter profiles being specific for different cortical layers.Here we propose a segmentation approach based on the observation that layers often cannot be distinguished using intensity alone in many histological stains.For example, layers III and V are both darkly stained by SMI-32 due to the presence of pyramidal cell soma [13] and consequently cannot be reliably distinguished using staining intensity alone ( Fig. 1 ).Hence combining staining information with a biologically and geometrically consistent depth map (i.e.depth profile) from readily available tools may provide a convenient way to create laminar segmentation from histological data.

Image acquisition
Feline experiments were approved by the local state authorities and were performed in compliance with the guidelines of the European Community European Union legislation for the care and use of laboratory animals (EU Directive 2010/63/EU), ARRIVE guidelines and the German Animal Welfare Act (TierSchG).Coronal sections of feline brain (right hemispheres) comprising all auditory cortical areas were collected.SMI-32 staining was used to identify cortical pyramidal cells [14] ( Fig. 2 a).Sixty SMI-32 stained slices, each 50 μm thick were obtained at 250 μm intervals.Overlapping image tiles were acquired under light microscopy at ∼1 μm isotropic resolution.They were then merged into whole slices using TeraStitcher [15] .

Volumetric reconstruction
The volumetric reconstruction of a curved object from two-dimensional slices is ambiguous without a shape prior [16] .Thus, this was done with the help of a volumetric feline brain atlas [17] consisting of a T1-weighted MR volume with 71 labeled brain regions.Histological images from four cats, two with normal hearing and two with congenital hearing loss, were volumetrically reconstructed in the following way.For each brain, the atlas was cropped to the brain region covered by the histology.Next the atlas was registered to the grayscale histology under a mean square error metric.Rigid histological slices were then aligned under a cost function, which simultaneously considered the error between adjacent slices and the MR volume [18] ( Fig. 3 ).
A region of interest consisting of the dorsal zone (DZ), primary (A1) and secondary auditory cortices (A2) was then isolated using the atlas labels followed by manual correction ( Fig. 2 b).This region was converted into a closed surface using restricted Delaunay triangulation [19] .This surface was cut along geodesic paths [20] into inner and outer surfaces at the gray-white and pial cortical surface respectively ( Fig. 2 c).

Equivolumetric depth
It has been hypothesized that cortical layer morphology varies due to equivolumetric deformation within cortical columns [ 21 , 22 ].In gyral columns, a constant volume ratio is maintained because outer regions are wider due to layers being thinner while inner regions are narrower due to layers being thicker.The opposite effect is observed in sulcal columns ( Fig. 1 ).Under this assumption, several methods for computing equivolumetric depth from images have been developed.These techniques can be broadly divided into two categories: voxel-based or surface-based.In voxel-based approaches, equivolumetric depth is approximated directly within the voxels of the image [23][24][25] while surface-based methods operate on triangulated meshes reconstructed from image data [26][27][28] .Here, we use a surface-based approach in which diffeomorphic registration of inner surfaces to pial surfaces was applied to create one-toone non-intersecting trajectories within the cortex [29] .These trajectories were defined on a time-domain  ∈ [0 , 1 ] from inner to outer surfaces with normality constraints enforced at each end.Equivolumetric depth was obtained through reparameterization of intermediate deformed surfaces to maintain constant volume along trajectories [30] .A corresponding depth image with range from 0 at the pia surface to 1 at gray-white matter boundary was then created by interpolating registration trajectories on a regular grid.Our goal is to approximate the parameters   of an -component Gaussian Mixture Model (GMM)    ∶  → ℝ that could have generated the joint histogram.It is defined as where    ∶  → ℝ is the  th zero-centered multivariate Gaussian probability density function with covariance   ∈ ℝ × .Furthermore   ∈ ℝ  are the means and   ∈ [0 , 1 ] are the weights such that ∑   =2   = 1 .Thus, our full set of parameters is defined as   = {  1 , … ,   ,  1 , … ,   ,  1 , … ,   } .

Algorithm
Given constant bandwidth, denoted by diagonal matrix ℎ = diag (ℎ 1 , … , ℎ  ) , the Filtered Kernel Estimate (FKE) is defined as where   (  ) =      ( −   ) −1   (  ) are filter functions and * denotes the convolution.The Alternating Kernel Method (AKM) adds components one at a time estimating new parameters θ which minimize the distance between the current GMM (1) and previous FKE (2).Thus, starting at  = 1 we find subsequent components by θ = argmin using gradient descent assuming that  θ0 =  (more details found in supplementary material).This process continues until the final number of components  =  max is attained.An example GMM and corresponding FKE after fitting to a joint histogram is depicted in Fig. 4 .

Initialization
Gradient descent optimization of Eq. ( 3) starts at an initial guess of parameters   = {  ,   ,   } ∪ { θ−1 } .These initial values were found automatically by modifying the "smart start " approach described previously [31] .First a new component with center   is added where the error  =  θ−1 −   −1 is greatest.In the original implementation, the histogram domain  was onedimensional and could be trivially divided into disjoint regions in which  > 0 .A new component was then inserted in the region of maximum error.Here in our multidimensional setting,  has a more complex topology and thus a different approach was needed.Let  ∶  → (ℝ ) be a centered-box kernel with unit sum.The convolution ( *  )(  ) can be interpreted as the average error within a region centered on  .Thus, the initial guess for the center of the new component is   = argmax  ∈ ( *  )(  ) .Once an initial   has been found, we assign the component a weight of   = ( *  )(  ) to ensure that the error is canceled out.Elements of covariance matrix   are then initialized to where (  )  ∈ ℝ and (  )  ∈ ℝ are the  th and  th elements of mean vector   respectively.

Segmentation
After convergence, a segmentation label image  ∶ Ω → {1 , . ..,  max } can be generated for multi-channel image .In the general case this can be done by But in specific case of laminar segmentation, this formulation can lead to unnatural results as subsequent labels should be contiguous (e.g.label 4 should always occur between labels 3 and 5).Suppose our components have been sorted by increasing depth from  = 1 to  max .Furthermore, suppose for a specific voxel  0 ∈ Ω we have  3 (( 0 ) ) = 0 .51 ,  5 (( 0 ) ) = 0 .49 and   (( 0 ) ) = 0 for  ∉ {3 , 5 } .Ideally a label of 4 should be assigned to this voxel since it is clearly between components 3 and 5.But instead the above equation would give us  ( 0 ) = 3 .Thus, we adopt where the ⌊⋅⌉ operation rounds values to the nearest integer.In the aforementioned example this would give us the desired result of  ( 0 ) = 4 .

Method validation
Grayscale histology from one of the deaf cats is shown in Fig. 5 .The first column depicts an axial slice through the region of interest followed by a coronal slice through the anterior ectosylvian sucal (AES) bank and then gyral crown.The second column shows the corresponding overlay of equivolumetric depth.The resulting multi-component segmentation is shown in the third column of ( Fig. 5 ).The 8th and 9th components clearly represented layer IV as shown in the fourth column, while an overlay of the manual segmentation is found in the fifth column.
Dice overlap [32] between manual and AKM segmentations were computed for all four brains.Average and standard deviation in these values across all four cats was 0.54 ± 0.03.Higher agreement 0.61 ± 0.03 was observed in gyral crown regions when compared to 0.50 ± 0.04 and 0.45 ± 0.10 found in the banks of the anterior (AES) and posterior ectosylvian sulci (PES) respectively ( Supplementary Table 1 ).Dice scores of < 0.7 were likely caused by well-known limitations of the Dice metric for smaller regions of interest [33] .Additionally, in sulcal banks, partial volume effects made segmentation more difficult due to lower resolution in slicing direction.Coronal slice planes were 0.25 mm apart while layer IV was typically 0.2 -0.3 mm thick.
Layer III/IV and IV/V interface surfaces were generated using restricted Delaunay triangulation for both the manual and AKM segmentations ( Supplementary Figs. 6 and 7 ).They were then compared using the Hausdorff distance.Since this metric is known for being sensitive to outliers [34] , particularly along the edge of the surfaces, the metric at 50th, 75th and 95th distance percentiles were calculated in additional to the traditional 100th percentile for the maximum distance.Across all brains, the mean and standard deviation for the 50th and 75th percentiles were 0.204 ± 0.043 mm and 0.315 ± 0.092 mm respectively ( Supplementary Table 2 ) as compared to a slice separation of 0.25 mm which suggests that layer IV can be reliably reconstructed in 3D.
A limitation of our segmentation approach is the reliance on an adequate spatial model in the form of an equivolumetric depth map.While in laminar segmentation a spatial model can be defined based on equivolumetric theory, this is not feasible in other biomedical segmentation problems.Furthermore, the class of kernel, Gaussian in our approach, is assumed a priori and limits the number of features that could be evaluated.Therefore, generalizing our method would require much more sophisticated methods with generalized kernels such as a U-Net trained on multichannel input images [35] .Qualitative agreement of algorithmic and manual segmentations demonstrates the potential for this technique to separating supragranular layers I-III from infragranular layers V-VI.Its utility in isolating individual cortical layers is more limited due to lower metric values.

Fig. 1 .
Fig. 1.Schematic of curved laminar coordinate system overlaid on SMI-32 histological sections.Roman numerals identify layer locations.Orange lines denote curved paths along cortical columns.Black and green lines denote gray-white and layer III-IV boundaries respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .Fig. 3 .
Fig. 2. (a) Original SMI Slice from a feline brain; (b) Region of interest consisting of DZ and A1 and A2 auditory cortical areas (marked in red); (c) Inner and outer cortical surfaces; (d) Cortical thickness (shown on the outer surface) are the lengths of the normal curves (shown in white) generated by constrained diffeomorphic evolution of the inner surface towards the outer one.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Results for a cat with hearing loss.Rows are axial slice, coronal slice through AES bank and coronal slice through gyral crown respectively.Columns are grayscale image, equivolumetric depth, MAKM 12-component segmentation and overlay on manual segmentation.Dashed lines indicate where axial and coronal slices intersect.